############
# packages #
############
library(PanelMatch)
library(plm)
library(dplyr)
library(ggplot2)
library(patchwork)
library(labelled)

###############
# set-up data #
###############
# read dataset
dat.anlsys <- read.csv("Party data for analysis.csv")
dat.anlsys <- as.data.frame(dat.anlsys)
# turn variables from integer to numeric
dat.anlsys <- dat.anlsys %>% mutate_if(is.integer, as.numeric)
# keep time and unit id as integer
dat.anlsys$period <- as.integer(dat.anlsys$period)
dat.anlsys$v2paid <- as.integer(dat.anlsys$v2paid)


##################################################
# roburstnes: 2 periods after change (Figure A4) #
##################################################
### matching and examination of set sizes
## moderate inclusiveness measure
# generate panel data for PanelMatch
dat.anlsys.dlgts <- PanelData(panel.data = dat.anlsys, time.id = "period",
                              unit.id = "v2paid",outcome = "v2pavote",
                              treatment = "prmrs.dlgts")
# perform matching 
PM.results.dlgts.2p <- PanelMatch(panel.data = dat.anlsys.dlgts, lag = 2, match.missing = TRUE, 
                               refinement.method = "ps.weight",
                               covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                 v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                 I(lag(gov.mmbr,1)) + I(lag(prmrs.dlgts.swtch.other,1)) + 
                                 e_gdppc + v2paallian, 
                               qoi = "att" , lead = 0:1, forbid.treatment.reversal = FALSE)
# examine set size
msets.dlgts.2p <- PM.results.dlgts.2p$att    
summary(msets.dlgts.2p)

## strict inclusiveness measure
# generate panel data for PanelMatch
dat.anlsys.mmbrs <- PanelData(panel.data = dat.anlsys, time.id = "period",
                              unit.id = "v2paid",outcome = "v2pavote",
                              treatment = "prmrs.mmbrs")
# perform matching 
PM.results.mmbrs.2p <- PanelMatch(panel.data = dat.anlsys.mmbrs, lag = 2, match.missing = TRUE, 
                               refinement.method = "ps.weight",
                               covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                 v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                 I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                                 e_gdppc + v2paallian, 
                               qoi = "att" , lead = 0:1, forbid.treatment.reversal = FALSE)
msets.mmbrs.2p <- PM.results.mmbrs.2p$att    
summary(msets.mmbrs.2p)

## set seed for boostrap se reproducability
set.seed(123)
## calculate estimates
PE.results.dlgts.90.2p <- PanelEstimate(sets = PM.results.dlgts.2p, panel.data = dat.anlsys.dlgts, confidence.level = 0.90)
PE.results.dlgts.95.2p <- PanelEstimate(sets = PM.results.dlgts.2p, panel.data = dat.anlsys.dlgts, confidence.level = 0.95)
PE.results.mmbrs.90.2p <- PanelEstimate(sets = PM.results.mmbrs.2p, panel.data = dat.anlsys.mmbrs, confidence.level = 0.90)
PE.results.mmbrs.95.2p <- PanelEstimate(sets = PM.results.mmbrs.2p, panel.data = dat.anlsys.mmbrs, confidence.level = 0.95)
## present numeric outcomes
summary(PE.results.dlgts.90.2p, verbose = FALSE)
summary(PE.results.dlgts.95.2p, verbose = FALSE)
summary(PE.results.mmbrs.90.2p, verbose = FALSE)
summary(PE.results.mmbrs.95.2p, verbose = FALSE)
## coefficient plots for both measures
# extract data frame from estimate
out.dlgts.95.2p <- summary(PE.results.dlgts.95.2p, verbose = FALSE)
out.dlgts.90.2p <- summary(PE.results.dlgts.90.2p, verbose = FALSE)
out.mmbrs.95.2p <- summary(PE.results.mmbrs.95.2p, verbose = FALSE)
out.mmbrs.90.2p <- summary(PE.results.mmbrs.90.2p, verbose = FALSE)
model.2p <- c("dlgts","dlgts","dlgts","dlgts",
              "mmbrs","mmbrs","mmbrs","mmbrs")
att.2p <- c(out.dlgts.95.2p[1], out.dlgts.95.2p[2], out.dlgts.90.2p[1], out.dlgts.90.2p[2], 
            out.mmbrs.95.2p[1], out.mmbrs.95.2p[2], out.mmbrs.95.2p[1], out.mmbrs.95.2p[2])
ci.upr.2p <- c(out.dlgts.95.2p[5], out.dlgts.95.2p[6], out.dlgts.90.2p[5], out.dlgts.90.2p[6], 
               out.mmbrs.95.2p[5], out.mmbrs.95.2p[6], out.mmbrs.90.2p[5], out.mmbrs.90.2p[6])
ci.lwr.2p <- c(out.dlgts.95.2p[7], out.dlgts.95.2p[8], out.dlgts.90.2p[7], out.dlgts.90.2p[8], 
               out.mmbrs.95.2p[7], out.mmbrs.95.2p[8], out.mmbrs.90.2p[7], out.mmbrs.90.2p[8])
ci.lvl.2p <- c(95,95,90,90,95,95,90,90)
cycle.2p <- c("t+0","t+1","t+0","t+1","t+0","t+1","t+0","t+1")
dat.PEresults.plot.2p <- data.frame(att.2p,model.2p,ci.upr.2p,ci.lwr.2p,ci.lvl.2p,cycle.2p)
dat.PEresults.plot.2p$cycle.2p <- factor(dat.PEresults.plot.2p$cycle.2p, levels = c("t+0","t+1"))

## plot results
PEresults.plot.2p <- ggplot(dat.PEresults.plot.2p, aes(y=att.2p,x=cycle.2p,group=cycle.2p)) +
  geom_point(size=2.5) +
  geom_linerange(dat=filter(dat.PEresults.plot.2p, ci.lvl.2p==95),aes(ymin=ci.lwr.2p,ymax=ci.upr.2p),linewidth=0.5) +
  geom_linerange(dat=filter(dat.PEresults.plot.2p, ci.lvl.2p==90),aes(ymin=ci.lwr.2p,ymax=ci.upr.2p),linewidth=1) +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_y_continuous(limits=c(-18,18), breaks=seq(-15,15,5)) +
  labs(y="Treatment Effect", x="Election After Treatment") +
  theme_bw() +
  facet_grid(. ~ model.2p, labeller = labeller(.cols = c(dlgts = "Moderate", mmbrs = "Strict"))) +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title = element_text(size=12),
    axis.text = element_text(size=11),
    strip.text.x = element_text(size = 12),
  )
## save plot (Figure A4)
ggsave(file="FigureA4.pdf", plot=PEresults.plot.2p, width=6, height=3.5)


####################################################################################
# robustness: matching by 3 lags in the DV and treatment (left panel on Figure A5) #
####################################################################################
### matching and examination of set sizes
## moderate inclusiveness measure
# perform matching 
PM.results.dlgts.3l <- PanelMatch(panel.data = dat.anlsys.dlgts, lag = 3, match.missing = TRUE, 
                                  refinement.method = "ps.weight",
                                  covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                    v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                    I(lag(gov.mmbr,1)) + I(lag(prmrs.dlgts.swtch.other,1)) + 
                                    e_gdppc + v2paallian, 
                                  qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# examine set size
msets.dlgts.3l <- PM.results.dlgts.3l$att    
summary(msets.dlgts.3l)

## strict inclusiveness measure
# perform matching 
PM.results.mmbrs.3l <- PanelMatch(panel.data = dat.anlsys.mmbrs, lag = 3, match.missing = TRUE, 
                                  refinement.method = "ps.weight",
                                  covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                    v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                    I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                                    e_gdppc + v2paallian, 
                                  qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# examine set size
msets.mmbrs.3l <- PM.results.mmbrs.3l$att    
summary(msets.mmbrs.3l)

## set seed for boostrap se reproducability
set.seed(123)
## calculate estimates
PE.results.dlgts.90.3l <- PanelEstimate(sets = PM.results.dlgts.3l, panel.data = dat.anlsys.dlgts, confidence.level = 0.90)
PE.results.dlgts.95.3l <- PanelEstimate(sets = PM.results.dlgts.3l, panel.data = dat.anlsys.dlgts, confidence.level = 0.95)
PE.results.mmbrs.90.3l <- PanelEstimate(sets = PM.results.mmbrs.3l, panel.data = dat.anlsys.mmbrs, confidence.level = 0.90)
PE.results.mmbrs.95.3l <- PanelEstimate(sets = PM.results.mmbrs.3l, panel.data = dat.anlsys.mmbrs, confidence.level = 0.95)
## present numeric outcomes
summary(PE.results.dlgts.90.3l, verbose = FALSE)
summary(PE.results.dlgts.95.3l, verbose = FALSE)
summary(PE.results.mmbrs.90.3l, verbose = FALSE)
summary(PE.results.mmbrs.95.3l, verbose = FALSE)
## coefficient plots for both measures
# extract data frame from estimate
out.dlgts.95.3l <- summary(PE.results.dlgts.95.3l, verbose = FALSE)
out.dlgts.90.3l <- summary(PE.results.dlgts.90.3l, verbose = FALSE)
out.mmbrs.95.3l <- summary(PE.results.mmbrs.95.3l, verbose = FALSE)
out.mmbrs.90.3l <- summary(PE.results.mmbrs.90.3l, verbose = FALSE)
model.l3 <- c("dlgts","dlgts","mmbrs","mmbrs")
att.l3 <- c(out.dlgts.95.3l[1], out.dlgts.90.3l[1], out.mmbrs.95.3l[1], out.mmbrs.95.3l[1])
ci.upr.l3 <- c(out.dlgts.95.3l[3], out.dlgts.90.3l[3], out.mmbrs.95.3l[3], out.mmbrs.90.3l[3])
ci.lwr.l3 <- c(out.dlgts.95.3l[4], out.dlgts.90.3l[4], out.mmbrs.95.3l[4], out.mmbrs.90.3l[4])
ci.lvl.l3 <- c(95,90,95,90)
dat.PEresults.plot.l3 <- data.frame(att.l3,model.l3,ci.upr.l3,ci.lwr.l3,ci.lvl.l3)
## plot results
PEresults.plot.l3 <- ggplot(dat.PEresults.plot.l3, aes(y=att.l3,x=model.l3,group=model.l3)) +
  geom_point(size=2.5) +
  geom_linerange(dat=filter(dat.PEresults.plot.l3,ci.lvl.l3==95),aes(ymin=ci.lwr.l3,ymax=ci.upr.l3),linewidth=0.5) +
  geom_linerange(dat=filter(dat.PEresults.plot.l3,ci.lvl.l3==90),aes(ymin=ci.lwr.l3,ymax=ci.upr.l3),linewidth=1) +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_y_continuous(limits=c(-13,13.5), breaks=seq(-10,10,5)) +
  scale_x_discrete(labels=c("Moderate", "Strict")) +
  labs(y="Treatment Effect", x="Inclusiveness Measure", title = "Match on 3-Lag DV") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title = element_text(size=12),
    axis.text = element_text(size=11),
  )

########################################################################
# roburstnes: including decade in matching (middle panel on Figure A5) #
########################################################################
### matching and examination of set sizes
## moderate inclusiveness measure
# perform matching 
PM.results.dlgts.decade <- PanelMatch(panel.data = dat.anlsys.dlgts, lag = 2, match.missing = TRUE, 
                                  refinement.method = "ps.weight",
                                  covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                    v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                    I(lag(gov.mmbr,1)) + I(lag(prmrs.dlgts.swtch.other,1)) + 
                                    e_gdppc + v2paallian + decade, 
                                  qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# examine set size
msets.dlgts.decade <- PM.results.dlgts.decade$att    
summary(msets.dlgts.decade)
## strict inclusiveness measure
# perform matching 
PM.results.mmbrs.decade <- PanelMatch(panel.data = dat.anlsys.mmbrs, lag = 2, match.missing = TRUE, 
                                      refinement.method = "ps.weight",
                                      covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                        v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                        I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                                        e_gdppc + v2paallian + decade, 
                                      qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# examine set size
msets.mmbrs.decade <- PM.results.mmbrs.decade$att    
summary(msets.mmbrs.decade)

## set seed for boostrap se reproducability
set.seed(456)
## calculate estimates
PE.results.dlgts.90.decade <- PanelEstimate(sets = PM.results.dlgts.decade, panel.data = dat.anlsys.dlgts, confidence.level = 0.90)
PE.results.dlgts.95.decade <- PanelEstimate(sets = PM.results.dlgts.decade, panel.data = dat.anlsys.dlgts, confidence.level = 0.95)
PE.results.mmbrs.90.decade <- PanelEstimate(sets = PM.results.mmbrs.decade, panel.data = dat.anlsys.mmbrs, confidence.level = 0.90)
PE.results.mmbrs.95.decade <- PanelEstimate(sets = PM.results.mmbrs.decade, panel.data = dat.anlsys.mmbrs, confidence.level = 0.95)
## present numeric outcomes
summary(PE.results.dlgts.90.decade, verbose = FALSE)
summary(PE.results.dlgts.95.decade, verbose = FALSE)
summary(PE.results.mmbrs.90.decade, verbose = FALSE)
summary(PE.results.mmbrs.95.decade, verbose = FALSE)
## coefficient plots for both measures
# extract data frame from estimate
out.dlgts.95.decade <- summary(PE.results.dlgts.95.decade, verbose = FALSE)
out.dlgts.90.decade <- summary(PE.results.dlgts.90.decade, verbose = FALSE)
out.mmbrs.95.decade <- summary(PE.results.mmbrs.95.decade, verbose = FALSE)
out.mmbrs.90.decade <- summary(PE.results.mmbrs.90.decade, verbose = FALSE)
model.decade <- c("dlgts","dlgts","mmbrs","mmbrs")
att.decade <- c(out.dlgts.95.decade[1], out.dlgts.90.decade[1], out.mmbrs.95.decade[1], out.mmbrs.95.decade[1])
ci.upr.decade <- c(out.dlgts.95.decade[3], out.dlgts.90.decade[3], out.mmbrs.95.decade[3], out.mmbrs.90.decade[3])
ci.lwr.decade <- c(out.dlgts.95.decade[4], out.dlgts.90.decade[4], out.mmbrs.95.decade[4], out.mmbrs.90.decade[4])
ci.lvl.decade <- c(95,90,95,90)
dat.PEresults.plot.decade <- data.frame(att.decade,model.decade,ci.upr.decade,ci.lwr.decade,ci.lvl.decade)
## plot results
PEresults.plot.decade <- ggplot(dat.PEresults.plot.decade, aes(y=att.decade,x=model.decade,group=model.decade)) +
  geom_point(size=2.5) +
  geom_linerange(dat=filter(dat.PEresults.plot.decade,ci.lvl.decade==95),aes(ymin=ci.lwr.decade,ymax=ci.upr.decade),linewidth=0.5) +
  geom_linerange(dat=filter(dat.PEresults.plot.decade,ci.lvl.decade==90),aes(ymin=ci.lwr.decade,ymax=ci.upr.decade),linewidth=1) +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_y_continuous(limits=c(-13,13.5), breaks=seq(-10,10,5)) +
  scale_x_discrete(labels=c("Moderate", "Strict")) +
  labs(y="Treatment Effect", x="Inclusiveness Measure", title = "Match on Decade") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title.x = element_text(size=12),
    axis.text = element_text(size=11),
  )

######################################################################
# roburstnes: 10,000 bootstrap iterations (right panel on Figure A5) #
######################################################################
### perform matching (identical to primary model)
## moderate inclusiveness measure
PM.results.dlgts <- PanelMatch(panel.data = dat.anlsys.dlgts, lag = 2, match.missing = TRUE, 
                               refinement.method = "ps.weight",
                               covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                 v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                 I(lag(gov.mmbr,1)) + I(lag(prmrs.dlgts.swtch.other,1)) + 
                                 e_gdppc + v2paallian, 
                               qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)

## strict inclusiveness measure
PM.results.mmbrs <- PanelMatch(panel.data = dat.anlsys.mmbrs, lag = 2, match.missing = TRUE, 
                               refinement.method = "ps.weight",
                               covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                 v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                 I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                                 e_gdppc + v2paallian, 
                               qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)

## set seed for boostrap se reproducability
set.seed(54321)
## calculate estimates (note that this will take longer to process than original model)
PE.results.dlgts.90.10k <- PanelEstimate(sets = PM.results.dlgts, panel.data = dat.anlsys.dlgts, 
                                         number.iterations = 10000, confidence.level = 0.90)
PE.results.dlgts.95.10k <- PanelEstimate(sets = PM.results.dlgts, panel.data = dat.anlsys.dlgts, 
                                         number.iterations = 10000, confidence.level = 0.95)
PE.results.mmbrs.90.10k <- PanelEstimate(sets = PM.results.mmbrs, panel.data = dat.anlsys.mmbrs, 
                                         number.iterations = 10000, confidence.level = 0.90)
PE.results.mmbrs.95.10k <- PanelEstimate(sets = PM.results.mmbrs, panel.data = dat.anlsys.mmbrs, 
                                         number.iterations = 10000, confidence.level = 0.95)
## present numeric outcomes
summary(PE.results.dlgts.90.10k, verbose = FALSE)
summary(PE.results.dlgts.95.10k, verbose = FALSE)
summary(PE.results.mmbrs.90.10k, verbose = FALSE)
summary(PE.results.mmbrs.95.10k, verbose = FALSE)
## coefficient plots for both measures
# extract data frame from estimate
out.dlgts.95.10k <- summary(PE.results.dlgts.95.10k, verbose = FALSE)
out.dlgts.90.10k <- summary(PE.results.dlgts.90.10k, verbose = FALSE)
out.mmbrs.95.10k <- summary(PE.results.mmbrs.95.10k, verbose = FALSE)
out.mmbrs.90.10k <- summary(PE.results.mmbrs.90.10k, verbose = FALSE)
model.10k <- c("dlgts","dlgts","mmbrs","mmbrs")
att.10k <- c(out.dlgts.95.10k[1], out.dlgts.90.10k[1], out.mmbrs.95.10k[1], out.mmbrs.90.10k[1])
ci.upr.10k <- c(out.dlgts.95.10k[3], out.dlgts.90.10k[3], out.mmbrs.95.10k[3], out.mmbrs.90.10k[3])
ci.lwr.10k <- c(out.dlgts.95.10k[4], out.dlgts.90.10k[4], out.mmbrs.95.10k[4], out.mmbrs.90.10k[4])
ci.lvl.10k <- c(95,90,95,90)
dat.PEresults.plot.10k <- data.frame(att.10k, model.10k ,ci.upr.10k, ci.lwr.10k, ci.lvl.10k)
## plot results
PEresults.plot.10k <- ggplot(dat.PEresults.plot.10k, aes(y=att.10k,x=model.10k,group=model.10k)) +
  geom_point(size=2.5) +
  geom_linerange(dat=filter(dat.PEresults.plot.10k, ci.lvl.10k==95),aes(ymin=ci.lwr.10k, ymax=ci.upr.10k),linewidth=0.5) +
  geom_linerange(dat=filter(dat.PEresults.plot.10k, ci.lvl.10k==90),aes(ymin=ci.lwr.10k, ymax=ci.upr.10k),linewidth=1) +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_y_continuous(limits=c(-13,13.5), breaks=seq(-10,10.5,5)) +
  scale_x_discrete(labels=c("Moderate", "Strict")) +
  labs(y="Treatment Effect", x="Inclusiveness Measure", title="10k Boostrap Iterations") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title = element_text(size=12),
    axis.text = element_text(size=11),
  )

## plot combined graph (Figure A5)
robust.combine.plot <- PEresults.plot.l3 + PEresults.plot.decade + PEresults.plot.10k +
  plot_layout(axis_titles = "collect")
ggsave(file="FigureA5.pdf", plot=robust.combine.plot, width=8, height=3.5)
